      Program plot_s
      !*****************************************************************
      include 'PPMA.Fort'
      Parameter(IMin=1,IMax=100,JMin=1,JMax=50)
      Dimension rho0(0:IMax,0:JMax),rho(IMin-1:IMax,JMin-1:4*JMax)
      Dimension r(IMin-1:IMax),teta(JMin-1:4*JMax)
      Dimension zp(-IMax:IMax,-IMax:IMax)
      Dimension zpt(-IMax:IMax,-IMax:IMax)

      pi=dacos(-one)
      do i=IMin-1,IMax
          r(i)=float(i)/float(IMax)
      enddo
      do j=JMin-1,JMax
          teta(j)=half*pi*float(j)/float(JMax) ! teta-symmetrical case
      enddo
      do j=JMin-1,4*JMax
          teta(j)=half*pi*float(j)/float(JMax) ! teta-symmetrical case
      enddo
      open(2,file='init_data')
      read(2,*) Re,crho,crhoe
      read(2,*) rho0
      crho0=one/rho0(0,0)
      do i=IMin-1,IMax
          do j=JMin-1,JMax
              rho(i,j)=crho0*rho0(i,j)
          enddo
      enddo

      do i=IMin,IMax
          do j=JMin,JMax
              rho(i,4*JMax-j)=rho(i,j)
          enddo
      enddo
      do i=IMin,IMax
          do j=JMin,JMax-1
              rho(i,2*JMax-j)=rho(i,j)
              rho(i,2*JMax+j)=rho(i,4*JMax-j)
          enddo
          rho(i,2*JMax)=rho(i,0)
      enddo

      dz=r(imax)/imax
      dx=r(imax)/imax
      dr=r(IMax)-r(IMax-1)
      dtheta=teta(JMax)-teta(JMax-1)
      do i0=-IMax,IMax
          do j0=-IMax,IMax
              z=i0*dz
              x=j0*dx
              r0=sqrt(z**2+x**2)
              if((i0.eq.0).and.(j0.eq.0)) then
              else
                  theta0=dacos(z/r0)
              endif
              if(x.lt.zero) theta0=2d0*pi-theta0
              xi1=r0/dr
              i1=xi1
              i1=min(IMax-1,i1)
              i2=i1+1
              xj1=theta0/dtheta
              j1=xj1
              j1=min(4*JMax-1,j1)
              j2=j1+1
              dr1=(r0-r(i1))/dr
              df1=(theta0-teta(j1))/dtheta
              zp(i0,j0)= (one-dr1)*(one-df1)*rho(i1,j1)
     +                  +(one-dr1)*df1*rho(i1,j2)
     +                  +dr1*(one-df1)*rho(i2,j1)
     +                  +dr1*df1*rho(i2,j2)
              zp(i0,j0)=log10(max(1d-10,zp(i0,j0)))
          enddo
      enddo
      do i0=-IMax,IMax
          do j0=-IMax,IMax
              zpt(i0,j0)=zp(j0,i0)   
          enddo
      enddo
      open(2,file='out.dat')
      write(2,*) 'x'
      write(2,*) 'z'
      write(2,*) -r(IMax),r(IMax),-r(IMax),r(IMax),zpt
      close(2)
      Call SM_Plot
      stop
      end

      Subroutine SM_Plot
      !*****************************************************************
      Parameter(IMin=1,IMax=100,JMin=1,JMax=50)
      parameter (IMax0=2*IMax+1,Jmax0=2*IMax+1,NCont=31)
      Dimension z(IMax0,JMax0),clev(NCont)
      character*90 ctime
      character*1 chx,chy

      open(2,file='out.dat')
      read(2,*) chx
      read(2,*) chy
      read(2,*) xmin,xmax,ymin,ymax,z
      xc=(xmin+xmax)/2
      yc=(ymin+ymax)/2
      zmin=z(1,1)
      zmax=z(1,1)
      do i=1,IMax0
          do j=1,JMax0
              zmin=min(zmin,z(i,j))
              zmax=max(zmax,z(i,j))
          enddo
      enddo
      zmax0=zmax
      zmin=-3.0
      do i=1,NCont
          clev(i)=zmin+(zmax-zmin)*(i-1.)/(NCont-1.)
      enddo
      call sm_EXPAND(1.4)
       call sm_device('postscript latypus')
c      call sm_device('xterm')
      call sm_graphics
      call sm_erase
      call sm_defvar('TeX_strings','1')
      call sm_limits(xmin,xmax,ymin,ymax)
      call sm_box(1,2,0,0)
      call sm_XLABEL(chx)
      call sm_YLABEL(chy)
      write(ctime,'(a8,a21,f5.2,a21,f5.2,a22,f5.2)') 
     +     'r_2=0.07',
     +     'log_{10}\rho_{min}=',zmin,
     +     'log_{10}\rho_{max}=',zmax, 
     +     '\Delta log_{10}\rho=',(zmax-zmin)/(NCont-1.)
      call sm_relocate(1.2*(xmin-xc),1.05*(ymax-yc))
      call sm_LABEL(ctime)
      call sm_EXPAND(1.0)
      call sm_defimage(z,xmin,xmax,ymin,ymax,IMax0,JMax0)
      call sm_levels(clev,NCont)
      call sm_contour
       call sm_hardcopy
      call sm_alpha
      return
      end
